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A numerical hydrodynamic model of a vibrated granular bed in 2D is elaborated based 
on a highly accurate Shock Capturing scheme applied to the compressible Navier- Stokes 
equations for granular flow. The hydrodynamic simulation of granular flows is numerically 
difficult, particularly in systems where dilute and dense regions occur at the same time 
and interact with each other. As an example of such a problematic system we investigate 
the formation of Faraday waves in a 2d system which is exposed to vertical vibration 
in the presence of gravity. The results of the CHD agree quantitatively well with event- 
driven Molecular Dynamics. 



1. Introduction 

Recently, an interest has developed towards the hydrodynamic simulation of complex 
granular fl ows, with the aim of understanding the details of transport in the contin- 
uum limit (|Goldhirschll2003 ) . Being both highly nonlinear and remarkably nonlocal, the 
mechanisms of granular transport in the dense limit are still obscure. The microscopic 
foundation of the hydrodynamic approaches of granular flow lies on the standard kinetic 
theory for molecular gases, except for the fact that dissipation is incorporated (via a 
constant or velocity dependent restitution) at the collision. While strictly, hydrodynamic 
descriptions cease to be valid near the close-packing limit, and at low restitutions, such 
models have however been used successfully in situations far from t heir supposed limits 
of validity, to desc ribe for inst ance shock waves in g r anular gases (iRericha et 
iBougie et a /] l2002h. clustering Jffijj fc Mazenkol 12003: iBrilliantov et alJl2QQ4h and coex- 



isting phases ([Meerson et q/.ll2003l 120041 ) . The difficulty of a hyd rodynamic description o f 
granular materials has been addressed in well reasoned terms in iGoldhirschl (fl999. 2001). 



One cannot pretend to overcome the problems of such a description, however, where it 
works exceptionally well, one should ask how, and why. 

In order to compensate for the energy lost in collisions, the typical mechanism of forc- 
ing, which is periodic vibration, induces not only the occurrence of various interesting 
phenomena but also the propagation of shock waves originating from the moving bound- 
ary into the system. The flow developed under these conditions is neatly supersonic and 
consequently, steep gradients arise in the hydrodynamic fields. For careful hydrodynamic 
simulations, it is then desirable to use shock-capturing methods to track the sharp fronts, 
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as well as high-order schemes which provide an accurate definition of the profiles. In the 
literature different methods have been applied to solve the Navier- Stokes equations for 
granular media, however, most of them lack the necessary analysis to discriminate the 
effects of the implicit or artificial diffusion added to handle supersonic flow. 

In order to show the possibilities of state-of-the-art hydrodynamic simulations (HD) 
applied to granular matter, in this paper we address a paradigmatic problem of granular 
dynamics: the observation of Faraday waves and the quantitative description/comparison 
of the instability via both hydrodynamic simulations based on WENO (Weighted Essen- 
tially Non- Oscillatory) schemes and the well established method of event-driven molecu- 
lar dynamics (MD). The difficulty that a hydrodynamic code faces when modeling such 
systems is twofold: the one introduced by the instability itself, regarding the use of proper 
parameters and the fine tuning of the code, and the one mentioned above, involving the 
propagation of shock waves across the granular bed under vibration and the discontinu- 
ities implied. 

The onset of patterns in vertically oscillated granular la yers has been nume rically 
investigated via particle and hydrodynamic simulations by iBougie et all ([2005). The 
numerical scheme employed there to solve the granular Navier- Stokes equations is based 
upon finite differences plus artificial diffusion. In the paper a thorough analysis is done on 
the role of fluctuations in the pattern inception and its dependence on the acceleration 
parameter, and the results mainly refer to the mean square displacement of the local 
centre of mass of the granular layer as an order parameter. Our analysis, being also 
focused on the comparison between particle and hydrodynamic simulations, is concerned 
with the whole set of hydrodynamic aspects related to the instability and how they are 
reproduced in order to assess the validity of granular hydrodynamic simulations in general 
in complex flow problems, which is not the object of study of the mentioned paper. 

The paper is organized as follows: we start in the next section summarizing the hydro- 
dynamic description of rapid granular flows based on kinetic theory. Section [3] is devoted 
to explain the details of the hydrodynamic code. In Section [H we describe the numerical 
hydrodynamic experiment and the molecular dynamics experiment; the results from both 
approaches are analyzed and compared in Section [5l In the conclusions we summarize 
our findings and expose a final comparison of both methods. 



2. Granular Kinetic Theory and Hydrodynamics 

As usual in granular kinetic theory, we follow the evolution of the number density of 
particles in phase space ( Brilliantov fc Poschelll2004l ). Given a 2D granular gas composed 



of small homogeneous disks of diameter a > and given (x, Vi), (x — an, V2) the states 
of two particles before the collision, where n G S 1 is the unit vector along the disks 
centers, the post-collision velocities are found by assuming that a fraction of the normal 
relative velocity is lost while its orthogonal component is unchanged. As a consequence, 
the post-collision velocities are obtained as 



1 V 

vi = -(v 1 + v 2 ) + T 

1 V' 

v^ = -(v 1 + v 2 )- T 



(2.1) 



where V = V - (1 + e) (V • n) n, V = Vi - V 2 and V = Vi - V 2 , with e being 
the coefficient of (constant) restitution. Let us denote by V* and V 2 the pre-collision 
velocities corresponding to Vi and V 2 . The Boltzmann equation for inelastic hard disks 
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gives the evolution of / (£, x, Vi) and it can be written as 

^ + (V 1 -V x )/ = ( x 2 Q(/,/) 



(2.2) 



where the collision operator is given by 
Q(f,f)=G(v)[ [ [(Vi-V 2 )-n] 

JR 2 JS 1 , 



-J/(vi)/(v;)-/(V!)/(v 2 ) 



dn d\2 



(2.3) 

in which dense gas effects are taken into account through the pair correlation function 
G(v), where v is the 2D volume fraction, i.e., v = jpcr 2 . The notation S]_ means that 
the above integral on n is taking over those values of n such that ((Vi — V2) • n) > 0. 
The moments of / allow us to compute the macroscopic quantities: the number density 



P (t, x 



)= / /(*,x, 



Vi) dV u 



the velocity field 



U(*,x) = - / Vi/(t,x,Vi)dVi, 
9 in 2 



and the granular temperature 

T(t,x) = ^ / |Vi-U(t,x)| 2 /(«,x,Vi)rfVi. 

Different formulas ar e proposed for G (y) in the literature but we choose the accurate 
formulas obtained bv lTorquatol (|l995h . More precisely, we use 



(2.4) 
(2.5) 
(2.6) 



G(v) 



1-0.436^ 



for ^ v < Vf 



(2.7) 



1 ~°- 43 % f for < v f < ^ c 

[l—UfY v c —v J 

where z// = 0.69 and ^ c = 0.82 are the freezing packing fraction and the random close 

packing fraction respectively^ 

By expansion methods ( Jenkins fc Richmanlll985a[ iGoldshtein fc Shapiro! 1 1995 ), the 
hydrodynamic equations are obtained for the 2D granular gas, 

dp 



dt 



+ V-(pU) =0 



p[ ^ + (U-V)U) =V-P + pF 
p(^ + (U-V)T)=-V-q + P:E- 7 



(2.8) 



representing the evolution of the number density of particles p(t,x), the velocity field 
U(t, x) = (Ui, U2) (£,x) and the granular temperature T(t, x). P = (-Pf?) is the pressure 
tensor, given in terms of the tensor E = (Eij) by 



with 



-p + (2/ii - /i 2 ) 

1 /at/,- <9t/,- 



(5y + 2//2-EV 



(2.9) 
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The principal pressure is related to the density and the temperature through the equation 
of state 

p = /9e[l + (l + e)G(i/)], (2.11) 
correction to the u s ual per f ect gases law, p = pe with e = T in 2D, pr o posed by 



correction to tne u s ual per l ect gases law, p = pe witn e = l m ZD, pr oposed by 
Jenkins fc Richmanl (|l985ah : iGoldshtein fc Shapirol (|l995h : iBrilliantov et all |2QQ4h to 
incorporate dense gas effects. The vector q = — xVT models the heat flux. The vis- 
cosities \±\ and /i2, the thermal conductivity x an d the cooling coe fficient 7 are density 
and tem perature dependent whose explicit expression we take from lJenkins fc Richman 
( 1985 oi l q). The bulk viscosity is given by 



V7T 



and the shear viscosity by 

^ = ^poTW 
the thermal conductivity is 

X = ^poT^ 
and the cooling coefficient is 



(1 - e 2 ) pT 3 ' 2 G{v). 



(2.12) 



(2.13) 



(2.14) 



(2.15) 



More invol ved hydro dynam ic models incorporate higher order terms in the gradients 
expansion dGoldhirs ch 2003), and more accurate expressions for the kinet ic coefficients 
which include high density corrections (iGarzo fc Duftv 1999t iLutskol 2005 ): in principle 
they can be easily implemented but it is more illustrative to keep the approach as simple 
as possible. Finally, the vector field F represents the external force acting on the system, 
for instance gravity modulated by the piston movement changing to the piston reference 
frame. The math ematical validity of t he hydrodynamic approximation in the Euler case 
was discussed in lBobvlev et aD (|2000h . 



3. Details of the hydrodynamic code 

In order to numerically simulate the granular Navier- Stokes equations (|2.8|h we write 
them as corrections to a compressible Euler- type system in conservation form. In fact, 
we consider the equivalent system 

g + V.(,U)=0 

^^ + V-[p(U®U)]+Vp= V-P + pF (3.1) 
dW 

— + V • \UW] = - V • q + P : E + U • (V • P) - 7 + p(U • F), 
where the total energy W is given by 

W = pT+^p\U\ 2 =pe+^p\U\ 2 . (3.2) 
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The system (j3.ip can be rewritten as 



with x = (#1, #2) and obvious definitions for f (u), g (u) and S (u) and with the vector 
u = (i/i, 1^2, ^3, 1/4) whose components are given by 

tzi = p 

^3 = PU 2 ( 3 ' 4 ) 

"4 = + l -p |U| 2 = pe + l -p (Ul + /7 2 2 ) . 

As a consequence, equations (j3.ip have the structure of a system of nonlinear conser- 
vation laws with sources. Local eigenvalues and both local left- and right-eigenvectors 
of the Jacobian matrices f ' (u) and g' (u) are explicitly computable and included in the 
appendix for the sake of the reader. Since local eigenvalues and eigenvectors need the 
derivative of the function G{y) appearing in the equation of state (|2.1ip . then G(v) is 
smoothed out around the freezing point volume fraction Vf. 

The numerical method consists in applying a high-order shock-capturing scheme for 
the Euler part of this system, i.e., the left-side of (|3.3p . while the terms in the right-hand 
side S (u) are approximated by second order central finite differences. 

The convective terms ^f^f (u) and ^f^-g (u) are approximated by a fift h-order finite 
difference characteristic- wise WENO method in a uniform grid following Uiang fc !Shu| 



(1996). This scheme is a well-known shock-capturing and high-order method for nonlinear 



conservation laws. Extensive benchmarks have shown this method to be particularly well 
adapted to co ntrol oscilla tions in shock regimes in classical Euler equations for gases, see 
the survey by|Shu| (|l998h . 

We first restrict ourselves to the one dimensional procedure, working in a dimension- 
by-dimension fashion for reconstructing each of the fluxes in Eq. (|3.3p . This means that 
when approximating (f (u)), for example, the other variable x 2 is fixed and the ap- 
proximation is performed along the x\ line. Given a uniform grid for any of the spatial 
variables denoted by r, let us call \ii(t) the numerical approximation to the point value 
u {ri , t) . The corresponding convective term is approximated as 

ilf( u ) ^ ~ 2 (3.5) 
dr v ' Ar K ' 

where f - , 1 is the numerical flux. 

2 

Let us briefly describe the general characteristic-wise finite difference procedure with 
flux splitting and flux reconstruction. At each fixed time t, we find the numerical fluxes 
f i+ 1 in the following way: 

Step 1. We first compute an average or intermediate state u i+ 1 using Roe's formula 
as inlShul (fl998h : iKamenetskv et all (|2000h . 



Step 2. We compute the matrices of right R and left R 1 eigenvectors and eigenvalues 
A/ ^u i+ i^, I — 1, . . .4, of the Jacobian matrix f ^u i+ i^ from the formulas postponed 
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to the appendix to obtain 

R = R(u, + i) 

R-^R- 1 ^) (3.6) 

A = A ( u i+^) = dia s ( u *+0 • 

Step 3. We transform the values and f (u^) in the potential stencil of the approx- 
imation of the flux to the local characteristic coordinates by using the left eigenvectors: 

Vi = R 1 ^ and h; = R _1 f (u*) . (3.7) 

Step 4. We use global Lax-Friedrichs flux splitting with the transport coefficient com- 
puted from the maximum, in absolute value, among the above computed eigenvalues to 
get from 

h± = l(h,±av,), (3.8) 



2 

where the dissipation parameter a is given by 



a = max 

u 



Step 5. We compute the decomposed fluxes h+ k and h , at the middle points using 
a high-order reconstruction explained below. 

Step 6. We transform back the computed fluxes to the physical space by using the 
right eigenvectors 

and finally, we add them to obtain the final numerical flux 

f i+i =f+ + i+fr +r (3-ii) 

Instead of Roe mean and Lax-Friedrichs flux splitting formulas, we can use more so- 
phisticated and improved methods avoiding numerical diffusion but in this case we do 
not need them since Navier-Stokes terms will have a regularizing effect anyhow on shock 
waves in the flow. 

In the WEN05 method, reconstructions of the fluxes in Step 5 are done using a non- 
linear nonlocal convex combination of three different approximations of the flux by three 
different local stencils using the upwinding of the flux. The contribution of each approxi- 
mation for the computation of the final value of the flux depend on the local smoothness 
of the function measured by certain smoothness indicators based on the local divided 
differences. If the function is smooth the resulting approximation is fifth order. The idea 
behind this weighted approximation is that those stencils including a discontinuity or 
high gradient will receive almost zero weight in the approximation and thus, oscillations 
will be avoided while keeping high-order approximation in smooth parts of the flow. 

More precisely, let us denote by h^ +1 j 2 l the l-ih component of the local numerical flux 

x+i/2 and by hf l the Z-th component of the local flux h^. Since the computations below 
are analogous component by component and fixed for the chosen upwinding, we will avoid 
the sub- and superindex "j for notational simplicity. Then, we obtain the numerical flux 

K+i/2 by 

K+l/2 = ^1^+1/2 + ^ill/2 + ^i+l/2 ( 3 ' 12 ) 
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where h^™\, 2 are the three third order fluxes on three different stencils given by 



^+ ) 1/2 = -^-i + ^ + ^+i, (3-13) 



1, 7, 11. 



•i+i/2 - 3^ g^+i e +2 ' 
and the nonlinear weights o; m are given by 



CJ m — 



Ez=i &i 

ll 



(3.14) 



with the linear weights ji given by 



1 3 3 

71 = ^, 72 = ^, 73 = ^, (3.15) 

and the smoothness indicators (3i given by 

13 9 1 9 

ft = 12 " 2hi ~ x + ^ + 4 " 4ft< " 1 + 3hi > 

P2 = § (fti-i - 2/2, + /\ z+1 ) 2 + i (^_i - /\,+i) 2 (3.16) 
& = § (ft< - 2/i i+ i + h i+2 f + i (3/i, - 4/^ +1 + /i i+2 ) 2 . 

Finally, e is a parameter to avoid the denominator to become and is taken as e = 10 -6 
in the computation of this paper. The computation for approximating h^_-jy 2 is obtained 
with the alternat e upwindin g by mirror symmetry with respect to the index We refer 
to the survey by Shu (1998) for further details about WENO reconstruction procedures, 
smoothness indicators, benchmarks and references for different applications. 

Finally, the resultin g ODE system is so lved in time by a 3rd-order TVD Runge-Kutta 
scheme introduced by Shu fc Osherl (|l998h imposing a local CFL condition at every time 
step, ensuring that the numerical speed is not less than the the largest eigenvalue of 
the local Jacobian matrix. This CFL condition becomes more restrictive as we approach 
the close-packing volume fraction. Although viscosity and thermal conductivity terms 
are treated as sources, they are taken into account in the CFL condition since diffusion 
effectively paces the time stepping of the code. This fact is due to the dependence of 
the viscosity and thermal conductivity on the volume fraction, resulting in diverging 
coefficients for both close-packing and vacuum limits. Both limits are avoided in the 
actual computation shown in this paper by establishing a limit from above and below of 
the volume fraction. These limits are 10 -4 and 99.99% of the maximal volume fraction. 



4. Numerical experiments 

4.1. Structures in vertically oscillated granular systems 

When granular matter is exposed to vertical oscillations in the presence of gravity, one 
observes spatio-temporal stru ctures at the free surface. Although this observation was 
reported as early as 1831 by iFaradavl the effect was studied only in recent ye a rs sys - 
tematically and a variety of patterns was observed, e.g. by EST & Behringer (1993); 





Figure 1. Periodic wave- like pattern in a vertically vibrated granular layer. The period of the 
pattern is twice the period of the excitation. The figure shows a two-dimensional MD simulation 
of 3000 particles at A = 4.02, / = 3.5 (details see below). Only a part of the system is shown. The 
dashed line indicates the lower amplitude of the oscillation. An animated sequence is provided 
as supplementary material. 



van Doom fc Behringer] (|l997h ; iMelo et all (|l994i Il995h ; [Umbanhowa r et all (Il996lll998h: 
Metcalfe et all ( 19971 ) and theoretically explained, e.g., bv lRothman ( 1998 ): Tsimring fc Aranson 



1997h . 

Here we focus on sub harmonic instabilitie s in a vertically vibrated layer of granular 
matter, first reported bv lDouadv et q/J ( 1989h and furth e r characterized in depe ndence on 
the parameters of the vibration by Clemen t et all ( 1996h : IWassgren et all (1996). Above a 
certain amplitude of acceleration, T = Auo 2 / g (A is the amplitude, uo the frequency of the 
sinusoidal excitation, g is gravity), a periodic wave-like pattern appears with frequency 
uj/2. Although the condition T > 1 is n ot a necessary condition for pattern formation 



( Poschel et al. 2000l : Renard et a/.|[200ll ). in all experimental observations the onset of 



the pattern formation was observed for T > 4. 

These waves could be re produced in quantit a tive agreement with e xperiments by nu- 
merica l MD sim ulations by Luding et al. (Il996l): iLudmd (Il997 ) and bv lAoki fc Akivama 
( 19961 ) (see also iBizon et all 19971 : Aoki fc Akivama 19971 ) and were also found from a 
linear stability analysis of an oscillating granular l ayer, modeled as an isothermal in- 
compressible fluid with vanishing surface tension ( Bizon et all 1 1999] ). Figure [T] shows 
snapshots of a Molecular Dynamics simulation. 

A computational fluid mechanics simulation of such a system represents a hard problem 
since the granular material occurs in dense state (e.g. snapshots 2 and 3 in Figure [1]) 
and in dilute state (e.g. snapshots 4 and 9 in Figure [Tj). A universal CFD computation 
should be capable to describe the system in both states dense and dilute. 

The described effect and its dependence on the system parameters was well studied in 
the literature; in the present paper it serves as an example showing that our CFD sim- 
ulation is in quantitative agreement with MD simulations and, thus, capable to describe 
time-dependent behavior of granular systems in dense and dilute state. For simplicity 
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and better visual representation, here we restrict ourselves to a two-dimensional system, 
the generalization to 3D is conceptually simple although computationally expensive. 

4.2. Setup of the Molecular Dynamics system 

Our reference system is an event-driven Molecular Dynamics simulation as shown in 
Figure [TJ In event-driven MD, each of the N particles of the system moves along a 
parabolic trajectory due to gravity g. This motion is interrupted by binary collisions, 
where each of the collision partners changes its velocity according to Eq. (|2.ip . Given 
particle positions and velocity at time t, the time of the next collision in the entire 
system can be computed analytically, therefore, in event-driven simulations, the time 
progresses in irregular steps, from event to event, i.e., from collision to the next collision. 
Unlike force-based MD, where Newton's equation of motion is solved numerically by 
evaluating the interaction forces, event-driven MD does not solve Newton's equation 
explicitly. Instead, Newton's equation is implied in the coefficient of restitution which 
relates the relative particle velocity after a collision to the velocity before the collision, 
Eq. (|2.ip . that is, the coefficient of restitution itself which, in general, is a function 
of material parameters and impact velocity, is a result of solving Ne wton's equation 
( Brilliantov et a/.lll996l: iSchw ager fc Poschell [l998: iRamfrez et !2[l999). In the present 



paper, we assume that the coefficient of restitution is a material constant and adopts the 
value e = 0.75. 

The main precondition for applying event-driven MD is the assumption of binary colli- 
sions which implies hard particles. Although conceptually simple, an efficient implemen- 
tation of event- driven MD algorithm s is by far more complex than force-based algorithms 
(see, e.g. JPoschel fc SchwagerlEoQ5[ ). To compare particle dynamics with hydrodynamics, 
it is necessary to employ event-driven MD since both MD and HD simulations are based 
on the same equations (j2.ip using both the concept of binary collisions and the coefficient 
of restitution. In general, there is no simple correspondence between force-based MD and 
hydrodynamics. 

The disadvantage of event-driven MD is its restriction to binary and thus, instan- 
taneous collisions. While for dilute granular gases this assumption is well justified, for 
systems with gravity, for mathema tical reasons we cannot e xclude tree-particle interac- 
tion, also called inelastic collapse (Mc Namara fc Youngll 19931 ). The re are numerical tricks 
to av oid or retard the inelastic collapse in MD simulations (e.g., iLuding fc McNamara! 
1998), provided density is small enough, however, for systems with gravity and reflecting 
oscillating bottom wall, at low frequency event-driven simulations become problematic: 
Consider a ball falling vertically from a certain height. At each collision with the wall 
its energy is reduced by the factor e 2 . Then, in finite time the ball comes to rest, that 
is, it stays in permanent contact with the wall which violates the condition of instan- 
taneous contacts. The same argument applies for slowly oscillating walls, therefore, our 
simulations fail for low oscillation frequency. 

To compare event-driven simulations with CFD results, we perform simulations of an 
assembly of N particles of diameter a = 10 mm, colliding due to Eq. (|2.ip with coefficient 
of restitution e = 0.75. The width of the system is L. After the transient time when the 
pattern is time-invariant (except for its period), we recorded the horizontal positions of 
all particles for 50 periods. The number of peaks N p of this pattern was then determined 
by visually inspecting the histogram of this data. In addition, we determined the wave 
lengths A of the pattern by Fourier- analysis, which agrees well with the values for N p . 

Table [T] shows the parameters of the performed event-driven MD simulations and the 
results for N p and A, extracted from Fourier analysis. The standard deviations dX account 
for the spread of the wavelength, which at very low amplitudes / accelerations is specially 
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Table 1. Performed MD simulations. N - number of particles, / = 2tiuj - frequency of the 
oscillating bottom wall (in Hz), A - amplitude of the oscillation, L - system width, H - number 
of particle layers. The observed patterns are characterized by N p - number of peaks, A - wave 
length, dX - standard deviation of the measured wavelength over 50 cycles. All lengths are given 
in units of the particle diameter. 
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high. The range of amplitudes A and accelerations T investigated is devised to capture 
the development of the instability; the width of the system L is chosen to ensure that 
a sufficient number of wavelengths fits the simulation window - and thus the errors in 
the determination of A are diminished. These results will be discussed below in Sec. [5] in 
comparison with the results of the computational hydrodynamics. 

Hydrodynamic fields such as density, velocity and temperature were generated from 50 
cycles of the MD simulation at amplitude A = 5.6 a (a: particle diameter) and frequency 
/ = 3.5 Hz, with the aim to compare them with the respective fields from the HD 
simulation. For this purpose we used a 150 a x 50 a grid. The positions and velocities 
of all particles were recorded for 50 periods at a rate of 50 fixed equidistant phases 
per period. In order to generate hydrodynamic fields, phase averaging was performed 
over particle positions and velocities with a spatial coarsegraining function (|Goldhirschl 
1999h (f)(r) oc (a/2 — |r — where 6 is the Heaviside step function and r^(t) the 

instantaneous location of the center of a particle. The cell size for averaging fits about 
3x1 particles, which leaves us with a typically mesoscopic procedure. 
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4.3. Setup for the hydrodynamic system 

We solve the full set of equations in 2D (|2.8p supplied with the kinetic coefficients given by 
Eqs. (|2.12p - (|2.15p and the equation of state (|2.1ip on a 150 x (50 up to 100) rectangular 
grid. The initial state is prepared by leaving the material settle down under the action of 
gravity to form a deposited layer of 6 grains. The depth of the granular bed is imposed by 
assuming a close random packing density in the deposited layer. As in MD simulations, 
the diameter of the grains is a = 10 mm and their restitution e = 0.75. 

The domain in which the hydrodynamic equations are solved varies in the range of 300 
to 600<j along the x\ direction, and from 30 to 100a along the vertical or x^ direction 
for the range of amplitudes A and frequencies / sampled in the simulations (1.5 to 15a 
in amplitude, 2 to 6 Hz in frequency). This way we make sure that the scale of the 
pattern is correctly modeled and a sufficient number of pattern wavelengths is captured. 
The choice of the vibration parameters is adjusted to cover the bifurcation region, that 
is, the region where the granular layer destabilizes from the homogeneous (along the x\ 
direction) solution. In general, we use exactly the same values for the parameters as in 
the MD simulations. 

We use periodic boundary conditions along the x\ direction, while at the top and the 
bottom we set reflecting boundaries. The top wall is set at a distance where reflections 
are of minor importance (the packing fraction at those heights is always of the order of 
10 -4 ). The hydrodynamic fields (density, velocity components, temperature) are stored 
for subsequent inspection without any additional treatment. Since the homogeneous so- 
lution tends to be very stable in the absence of random sources, we numerically perturb 
the layer density with a collection of modes (the first 20 with largest wavelengths), with 
random phases and amplitudes. The latter are selected with a maximum in the order of 
10 -3 times the value of the local packing fraction of the deposited material. This small 
perturbation is sufficient to quickly observe the growth of the m ode selected by the system 
which gives rise to the development of the Faraday waves (see [supplementary materiall ) . 



5. Results 

Unlike particle simulations, in hydrodynamic simulations the patterns tend to be very 
stable and regular for small as well as for moderate amplitudes. The growth of the 
sele cted wavelength develops within 10 cycles after the perturbation has been imposed 
(see supplementary materiall ) , and one observes only minor changes after this time. The 



periodicity of the pattern is twice the shaking period, as corresponds to the typical period 
doubling of this instability. Figure [2] shows the appearance of the pattern in a sequence of 
eight frames which covers two periods of shaking and the corresponding frames generated 
from the MD simulation. 

The main difference between the HD and MD sequences in Figure [2] is the different 
landing time of the granular bed. The fourth and the eighth snapshots are taken when 
the material visibly starts to deposit, which happens at different times in HD and MD 
simulations. This deserves a careful discussion. Indeed, one of the typical features ac- 
companying the instability is the gap formed between the material and the bottom wall, 
induced by an acceleration of the plate at T > 1. In real experiments as well as in MD 
simulations - and leaving aside the effects of particle noise, there is a true empty space 
opening and closing periodically. In HD simulations the density is never very close to 
zero at the bottom wall. This fact has to be attributed to the implemented boundary 
condition. More will be said in this respect later in this section. 

The bifurcation diagram is shown in Figure El from which one can easily observe the 
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Figure 2. Eight pictures of the density as obtained from the numerical solution of the hydrody- 
namic equations (HD, left) and the corresponding ones from the molecular dynamics simulation 
(MD, right). The times for the HD snapshots are: 0, 0.25T, 0.5T, 0.75T, T, 1.25T, 1.5T, 1.75T, 
with T being the period of the oscillation. The MD snapshots were taken at times: 0, 0.25T, 
0.5T, 0.86T, T, 1.25T, 1.5T, 1.86T. The initial phase is arbitrary, but common for both sets. For 
most of the period, the dynamics observed in HD and MD simulations is completely equivalent; 
the difference in the fourth and eighth phases is due to different landing times of the granular 
bed (see text). The MD figures are obtained after averaging particle positions over 50 cycles, 
while the density field on the left is the raw output. The frequency of vibration is / = 3.5 Hz and 
the amplitude A = 5.6 diameters. Complete sequences are available as supplementary material. 



extent of the agreement that both types of simulations can provide, with results which 
quantitatively diverge by no more than 20%. Close to the bifurcation point though, there 
is an important exception: According to the hydrodynamic simulations, the pattern starts 
to form between 1.55 and 1.64 diameters of amplitude, since below A = 1.64a it is not 
observed for any reduced acceleration T = A(27rf) 2 /g in the range from 2.0 to 2.75, where 
g is the acceleration of gravity. For the MD simulations, however, one obtains a few data 
points in the region of very low amplitudes giving nonzero wavelengths. Nevertheless 
the errors associated with these few MD measures are particularly large, coming from 
the procedure to determine the characteristic pattern wavelength in MD simulations: 
due to the general contribution of particle noise, the wavelength was extracted from the 
histogram of the Fourier analysis of the particle positions over 50 periods. Otherwise, in 
the low amplitude regime, no pattern would be guessed from the naked-eye inspection of 
the MD sequences. Therefore, the typical histograms contain a wide spread of wavelengths 
in this regime, which reflects in the standard deviations (dX in Table HJ plotted in Figure 
[3] as error bars. The duplicate data points in some of the hydrodynamic samples refer to 
situations where a child peak develops or disappears (due to finite size effects introduced 
by the periodic boundaries) and as a result, both wavelengths can be observed in the 
course of one simulation. This tends to happen in the large amplitude regime, where the 
pattern is less stable. 

In Figured! we show a 2-period sequence of the vertical profiles of the packing fraction 
v and the rescaled internal energy pT /erg, as a function of height in diameters. For 
the plots we have selected one of the locations in Figure [2] where a valley /cusp develops 
alternatively. The first snapshot (a) shows the granular material already being compressed 
by the incoming bottom wall. As one can appreciate, the extent of the layer (solid line) 
is small so the profile corresponds to the location of a valley at this time. The thermal 



Pattern Formation in Vertically Oscillated Granular Disks Layers 



13 



100 



80 



60 



40 



20 




o 


r = 


2.75 (MD) 




□ 


r = 


2.5 (MD) 




o 


r = 


2.25 (MD) 




A 


r = 


2 (MD) 




• 


r = 


2.75 (HD) 




■ 


r = 


2.5 (HD) 




♦ 


r = 


2.25 (HD) 




▲ 


r = 


2 (HD) 





5 10 15 

amplitude (in diameters) 



20 




J \ 2 3 4 5 
amplitude (in diameters) 



Figure 3. The wavelength of the pattern as a function of the amplitude of shaking, for different 
reduced accelerations Y. Open symbols: Numerical solution of the hydrodynamic equations; full 
symbols: Event-driven molecular dynamics simulations. On the right: A zoom of the figure at 
low amplitudes. Error bars are drawn only for MD data points. 



energy is high and has started to propagate upwards after the impact. In Figure ffl(b) 
the layer is thickening while at the peaks is thinning, as (f) shows. In (c) the energy 
has already propagated to the top, the material is suspended and rather fluidized, and 
the cusp is visible. In (d) another impact is taking place, the density is growing at the 
bottom, and one more shock wave has been generated whose energy extincts as the peak 
redissolves to feed the closest locations where new cusps are forming - indeed, in (e) 
the material is close packed at the bottom, and thus density can only further decrease 
by pressure gradients. From the comparison between (c) and (g), one can see that the 
densities are very different in peaks and valleys, and the shock progresses thus differently. 
The careful analysis of the shock wave pr opagating through a (homogeneous) vibrated 
bed has been done bv lBougie et all ( 20021 ) and here suffices to point at the fact that the 
flow is indeed supersonic, with Mach numbers in the range from to 10. See the Figure 
[5] to observe the typical oscillations of the Mach number at a valley / cusp. The highest 
values of the Mach number are achieved in the dilute phase but not far away from the 
clustered phase; the figure shows the Mach number recorded at the height where the 
packing fraction is 5%, as a convenient measure of the interface between the dense and 
rarefied material. The convective motion along the x\ direction which removes material 
from the dense towards the dilute locations will be analyzed later. 

For the MD simulation, the time-dependent thermal energy has been extracted by av- 
eraging over 50 cycles the corresponding phases of the local fluctuational kinetic energy, 
as defined by Eq. (|2.6p . Afterwards an average over equivalent locations has been per- 
formed. The profiles are shown in Figure [6j Comparison between Figs. [Hand [6] indicates 
that the quantitative difference in the temperature profiles is small. However, the hydro- 
dynamic profiles show a sharper shock wave, as compared with their corresponding MD 
counterparts. This feature is related to the residual heating observed in MD simulations 
at heights where the granular layer is sufficiently dilute (very vi sible for example in Figur e 
Et). Movies showing the complete sequences are available as (supplementary material! ) . 
From their inspection one can realize that the maximum values of the temperature do not 
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Figure 4. Hydrodynamic simulation: A 2-period sequence of the vertical profiles of the packing 
fraction v (solid line) and the rescaled internal energy pT/ag (dashed) as a function of height, 
in diameters, for the same times as the snapshots in Figure [2j left). 



differ very significantly, whereas the MD profiles appear systematically thicker. This phe- 
nomenon might be attributed in part to the mesoscopic averaging procedure employed 
here. It is certainly known, on t he other hand, th at hydrodynamic variables can be scale 
dependent in granular systems (|Goldhirschll2Q0lh due to the lack of scale separation. So 
even if the number of samples (cycles) could in principle be increased, and therefore the 
mean fields better defined, the local noise will persist at scales of the order of one particle 
(which is the grid size used here for averaging the MD sequences), particularly in dilute 
regions. The number of samples is practically limited by the fact that the MD pattern 
is not strictly stationary for very long times. It is not desirable either to increase the 
grid size for averaging, for obvious reasons: Here, where the pattern is of the size of a 
few grains, the fields obtained from MD simulations would lose structure and definition. 
Oppositely, the granular Navier-Stokes equations (|2.8|) do not account for such a source 
of fluctuations. However, and even with these restrictions, one can conclude that good 
quantitative agreement can be achieved with mesoscopic statistics: a manageable number 
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Figure 5. Solid line: Evolution of the Mach number over two periods for the frames in Figure 
[4] at a height (variable) where the packing fraction is 5%. Dashed line: height (in diameters 
divided by 10) corresponding to the values of the Mach shown. 



of samples (25) and a one-diameter characteristic cell size, along with a final average of 
means over corresponding locations. 

There is still another difference between the MD and HD profiles, and it is related to 
the gap formed when the granular layer takes off (Figure [2k and g). In MD simulations, 
during the stage of downward displacement of the wall, the entire layer is levitating, as 
one can observe from the corresponding frames in Fig. [6] by noting that the density at 
zero height is exactly zero. However, in HD simulations the density at the bottom is 
never very close to zero. It is particularly large at the base of the cusps, (Figure Ht), 
where the packing fraction is 0.19. Seemingly, the hydrodynamic layer is more sticky, or 
behaves more inelastically, than the MD granular bed. Besides that, in hydrodynamic 
simulations, the fact that there is some material stuck to the bottom during the take off 
anticipates the impact with the wall and the corresponding generation of the shock wave. 
This is the source of the mismatch between the landing times, to which Figure [2] refers. 

It is worth to analyze the velocity field and the convection pattern which accom- 
panies the Faraday waves. These are shown in values relative to the velocity of the 
plate in Figs. [7] an d [8] for hydrodynamic and MD simulations, respectively, and in the 



supplementary material In fact what is shown is the linear momentum field rather than 



the velocity, in order to focus onto the most relevant part of the system eliminating from 
the picture the almost empty regions. Comparing both panels we see again the effect that 
was pointed at above: in hydrodynamic simulations, the material never abandons contact 
with the vibrating plate. This feature has been observed befor e in hydrodynam ic simu- 
lations of vibrated granular layers without pattern formation (|Shattuckl 120051 ). Besides 
that, in Figs. [3 and [8] we observe the removal of material from the dense towards the loose 
regions after each impact of the wall, during the alternate formation of cusps. The scaling 
factor for the arrows - chosen identical for both figures, again reveal good quantitative 
agreement of the hydrodynamic simulations and MD. In HD we obtain slightly lower 
values of the linear momentum as compared with MD simulations, by a factor which is 
shown to be about 20% - the maximum value achieved by the linear momentum was 1.17 
in some units for the MD sequence; for the HD simulation was 0.935. We recall that sim- 
ple reflecting boundaries are used for the hydrodynamic simulations (with zero no-slip), 
and that particles are perfectly smooth in MD simulations. It is difficult to evaluate the 
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Figure 6. Molecular Dynamics simulation: A 2-period sequence of the vertical profiles of the 
packing fraction v (solid line) and the rescaled internal energy pT/ag (dashed) as a function of 
height, in diameters, for the same times as the snapshots in Figure fright) . 



effect of the choice of boundary conditions in more complex setups, but in view of the 
results, it seems clear that this undisturbed sliding along the plate looks like the only 
reasonable choice in this very simple situation. 



6. Conclusions 

Up to now, we have compared the results given by the two simulation methods con- 
sidered: The traditional approach by event-driven MD and the hydrodynamic approach 
by means of a CFD code. We have shown that qualitatively and quantitatively as well, 
adequately armored hydrodynamic codes can reproduce the features of particle simula- 
tions in complex flow problems dealing with transient strong shock waves and pattern 
formation, such as the Faraday instability here investigated. In particular, the appear- 
ance and wavelength of the pattern has been reproduced with a very good approximation 
as seen from the bifurcation diagram (Figure [3]). The hydrodynamic fields compare very 
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Figure 7. Hydrodynamic simulation: The linear momentum plotted over the density field in 
gray shades. The frames are taken at the same times as in Figure [2{left) and cover two cycles 
of the driving oscillation. The complete sequence is available as supplementary material. 
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Figure 8. Molecular Dynamics simulation: The linear momentum plotted over the density field 
in gray shades. The frames are taken at the same times in Figure [2j right) and cover two cycles 
of the driving oscillation. 
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well to the corresponding averaged particle fields (with the exception of the gap formed 
during the take off of the granular layer). The hydrodynamic temperature is shown to 
agree with its MD counterpart to which ensemble averaging at a typical grid size of one 
particle diameter has been applied. Last, the velocity field obtained in hydrodynamic 
simulations follows the periodic alternancy of the pattern and shares the type of struc- 
ture observed in the particle velocity field, if however giving somewhat smaller values. 
The discrepancies found along the analysis are of the order of 20%. 

It is worthwhile now to compare the efficiency of the two methods in general. In the 
first place, one does not reveal any secret saying that CFD codes are computationally 
expensive. On the other hand, having said that the scale of the gradients is one particle 
diameter, one has served the trouble: Effectively, the cell size cannot be greater than a 
few particles. At values of the Mach number of about 10 in some phases of the motion, 
inertia is the most prominent mechanism of transport. Diffusion operates at very small 
length scales, and is of less importance than inertia, but plays a role in some other 
phases of the motion, when the typical Mach number drops below 1 in the dense phase, 
and cannot be just neglected. In addition, explicit time schemes like the one used in 
this paper introduce a stability condition which turns more and more restrictive as the 
diffusion length scales (the size of the particles) become smaller. The result is that for 
submillimetric particles, most commonly used in real experiments, the same problem 
treated here turns unbearably stiff, as the typical timescale of the pattern does not 
decrease in the same correspondence. 

In any case one must resolve length scales of the typical size of one particle, and the 
ratio one cell w one particle is extremely unfavorable to computational hydrodynamics in 
terms of efficiency. Moreover, high resolution shock-capturing schemes like WENO add 
an extra cost in the computation which often amounts the largest portion of the total 
CPU time. There is some room for optimization though, in the sense that not everywhere, 
not all the time, such extraordinary effort to construct highly accurate numerical fluxes 
should be necessary, and thus the simulation time could this way be reduced. 

The question about boundary conditions still deserves an additional remark. Here the 
most simple choice has been made, with specular reflecting boundaries without any por- 
tion of no-slip in the tangential direction (the velocity component parallel to the plate 
is not disturbed by the boundary). Presumably, in order to mimic rough particles some 
partial slip would have to be introduced with parameters to be determined from MD 
simulations, although friction may trigger rotational mechanisms which are not contem- 
plated in the set of hydrodynamic equations considered here. 

Finally, one should ask why the agreement between molecular dynamics and hydrody- 
namics is here so good, so robust, when all ingredients are present to make it fail. At high 
densities collisions are not binary, fluxes are thus nonlocal and the expressions for the 
kinetic coefficients should no longer apply, not to mention the degree of inelasticity (the 
restitution is only e = 0.75) which invalidates the assumption of small gradients on which 
all kinetic theory expansions are based. It seems that there is only one reason which can 
explain the agreement: The dynamics is mostly Eulerian and details like constitutive 
relations and transport coefficients are not definitive. So much so that in practice, if one 
focuses the attention on modeling accurately the shock wave and introduces a sensible 
equation of state, the instability is "easily" reproduced, as in molecular dynamics, with 
the essential contribution of inelasticity. 
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Appendix A. Formulas for local variables 

Here, we summarize the formulas for the explicit eigenvalues and eigenvectors of the 
Jacobian matrix of the Euler part of (|3.ip in 2D, for general equations of state depending 
on the density and the enthalpy. As discussed in Section [3l the system (|3.ip can be 
rewritten as 



<9u d d 



(Al) 



with obvious definitions for f (u), g (u) and S (u) and with the vector u = (u±, u 2l ^3, U4) 
whose components are given by 
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where the pressure is assumed to be a function of density and enthalpy, 
The sound speed is given by 
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Given the auxiliary function is now H = (p + u±) ju\, the eigenvalues of the Jacobian 
matrix f / (u) are given by 
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and their corresponding right and left eigenvectors are 
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respectively, where b\ = p e / pc 2 s . For the fluxes associated with the x 2 -derivative, that 
is, 

u 2 u 3 



g x (u) = u 3 , g 2 (u) = 

2 ^ 

g 3 (u) = p H , # 4 (u) = — (p + u A ) ; 



(A 9) 



we just need to swap indices 2 and 3 in all formulas above, and the second and third 
components in all right and left eigenvectors. 
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